#!/usr/bin/perl
#use strict;
# use warnings;
use Data::Dumper;

#####################################################

my ($USAGE) =  "\n\n\t****** USAGE of panam.pl PROGRAM ******\n\n\n\n\tUSAGE: perl $0 <panam.ini file> <taxonomic_group> \n\n\n\n";
die "$USAGE" if ((scalar(@ARGV))< 1);
# my $option_file = $ARGV[0];
my $option_file = "panam.ini";
chomp($option_file);
my $taxo_group = $ARGV[0];
chomp($taxo_group);

die "\n\n\t=> Cannot find configuration file: $option_file.\n\n" unless (-e $option_file);
die "\n\n\t=> Configuration file, $option_file, appears to be empty!\n\n" if (-z $option_file);
open(OFILE, "<$option_file") || die "Cannot open input file : $option_file\n";

my $usearch;
my $NGS_id_Results;
my $panam_output = "panam_output";
# my $parsing;
my $query_seq; 
my $user_file;
my $preprocess_output = "preprocess_output";
my $seq_F;
my $seq_R;
my $primF;
my $primR;
my $trim =0;
#my $path = "bin/uclust3.0.617";
my $us_version = "3.0.617";
####
my $lib = "lib";
#my $fast = "$lib/FastTree";
#my $align = "$lib/hmmalign";
my $Reference = "Reference" ;
#my $build = "$lib/hmmbuild";
####
my $dom;
my %dom;
my $aff_norm;


while (<OFILE>) {
	chomp;
	my $line = $_;

# discard blank line
	if ($line =~ m/^\s*$/g) {
		next;
	}

# discard comment line
	elsif($line =~ m/^\s*#/g) {
		next;
	}
	
# get input options line
        else {
		
		if($line=~m/454_RUN_IDENTIFIER/gi) { 
			$line=~s/454_RUN_IDENTIFIER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$NGS_id_Results = $line; 
			}
	 		die "\n\tOutput directory for NGS analyses is not defined. Check panam.ini.\n\n" unless ( $NGS_id_Results ne "");
			unless (-d $NGS_id_Results) { mkdir $NGS_id_Results || die "\n\tCould not create Directory $NGS_id_Results $!\n"; }
		}

# 		if (-d $NGS_id_Results."/".$panam_output) { `rm -r $NGS_id_Results/$panam_output` } 
		if (!(-d $NGS_id_Results."/".$panam_output)) { mkdir $NGS_id_Results."/".$panam_output || die "\n\tCould not create Directory $NGS_id_Results/$panam_output $!\n"; }


		if($line=~m/QUERY_SEQUENCES/gi) {
			$line=~s/QUERY_SEQUENCES//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi; 
			if($line) {
				$query_seq = $line; 
			}
			chomp ($query_seq);
			# See if the query sequences file exist and contains something
			die "\n\n\tCannot find query sequence file: $query_seq. Check panam.ini.\n\n" unless (-e $query_seq);
			die "\n\n\tQuery sequence file, $query_seq, appears to be empty!\n\n" if (-z $query_seq);
			
			my $ligne1= `sed -n '1p' $query_seq`;
			my $ligne2= `sed -n '2p' $query_seq`;	
			if (($ligne1 !~ /^>/) or ($ligne2 !~ /^[atcguATCGU]/)) {
				die "\n\n\tYour file $query_seq does not seem to contain sequences in fasta format.\n\n";
			}		
		}

		if($line=~m/USER_FILE/gi) {
			$line=~s/USER_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$user_file = $line; 
			}
			die "\n\n\t User file value missed. Check panam.ini.\n\n" unless (($user_file ne ""));
			die "\n\n\t Invalid user file value!\n\n" unless (($user_file eq "own") or ($user_file eq "pooledSamples") or ($user_file eq "eachSample"));	
		}

		if($line=~m/FORWARD_PRIMER_NAME/gi) {
			$line=~s/FORWARD_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primF = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_NAME/gi) {
			$line=~s/REVERSE_PRIMER_NAME//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$primR = $line; 
			}
		}

		if($line=~m/FORWARD_PRIMER_SEQUENCE/gi) {
			$line=~s/FORWARD_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_F = $line; 
			}
		}
	
		if($line=~m/REVERSE_PRIMER_SEQUENCE/gi) {
			$line=~s/REVERSE_PRIMER_SEQUENCE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seq_R = $line;
			}
		}

		############### normalisation
		if($line=~m/NBR_SEQ_NORM/gi) {
			$line=~s/NBR_SEQ_NORM//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$aff_norm = $line; #print "aff --- $aff_norm\n";
			} 
		}
		
		#####################################

		if($line=~m/DOMAIN/gi) {
			$line=~s/DOMAIN//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) { 
				$dom = $line;			
				if (defined $dom ) { #push (@dom, $dom);
					$dom{$dom} = 1
				}
			}
			foreach my $d (keys %dom) {
				die "\n\n\t Domain name is not correct. Check panam.ini.\n\n" unless (($d eq "bacteria") or ($d eq "eukaryota") or ($d eq "archaea"));
			}	
		}
	}
}

die "\n\n\t Invalid number of sequences to pick for the normalization! Check panam.ini.\n\n" if ((defined $aff_norm ) and ($aff_norm !~ /^\d+$/));
if (!(defined $aff_norm)) {$aff_norm = '200' }


print "$Reference\n";

############################################# trimmed profiles path ######################################################

if (!(defined $primF) and (defined $primR))  {
	die "\n\n\t Primer forward not defined. Check panam.ini.\n\n";
}

if (defined $primF) {
	die "\n\n\t Primer reverse not defined. Check panam.ini.\n\n" unless (defined $primR);
	die "\n\n\t Primers' sequences missed. Check panam.ini.\n\n" unless ((defined $seq_F) and (defined $seq_R));
	die "\n\n\t Primers' sequences seem to be in the wrong format. Check panam.ini.\n\n" unless (($seq_F =~ /[ATCGUatcgu]/) and ($seq_R =~ /[ATCGUatcgu]/));
	$trim = 1
}

my @tax=();
my $Profiles;

if ($trim == 1) {
	#$Profiles = "Reference/Profiles_".$primF."_".$primR;
	$Profiles = "Reference/Trimmed_Profiles";
	if ((-d $Profiles)) {
		print "Trimmed profiles exist in $Profiles/\n";

	}
}

 
if ($trim == 0) {
	$Profiles= "Reference/Profiles";
}

# 
# ########################################################################################################
# 
# ########## hash reference sequences

#`mkdir $NGS_id_Results/$panam_output"/Similarity_annotation/Sorted_Sequences"`;


###################################

# ########## hash query sequences 
my %h; my $c;
open (F, $query_seq);
while (<F>){
	if ($_ =~ /^>(.*?)\n/) {
		$c=$1;
		$c =~ s/\s+$//; 
		
	} else {
		chomp ($c);
		$h{$c}.=$_ ;
	}
}
close F;


########### hash tax                            
my %taxeuk;
open (T, "Reference/Taxonomy") || die "can not open file";
while (<T>) {
	my @e = split (/\t/, $_);
	chomp($e[0]); chomp ($e[1]); chomp ($e[2]);
	$taxeuk{$e[0]}=$e[1]."\t".$e[2];
}
close T;

####################### sorting sequences

open (F, $NGS_id_Results."/".$panam_output."/Similarity_annotation/seq_results.uc") || die "can not open file" ;
my %hits;
while (<F>) {
	if (!($_=~ /\#/)) {
		my @t = split (/\t/, $_);
		if ($t[0] eq "H") {
			chomp ($t[8]) ,chomp($t[9]);
			$t[8] =~ s/\s+$//; $t[9] =~ s/\s+$//;
			$hits{$t[8]} .=$t[9]."/"; 

			if ($t[4] eq "-") {
				my $rev = reverse ($h{$t[8]});
				$rev =~ tr/ACGTU/TGCAA/ ;
				chomp ($rev);
				$h{$t[8]} =$rev;
			}
		}
	}
}
close F;

open (B, ">".$NGS_id_Results."/".$panam_output."/Similarity_annotation/best_hit_uc");

my %tax_hits; my %seq_hits;
foreach my $a (keys %hits) {
	
	my @tab = split (/\//, $hits{$a});
	my $i = 1;
	foreach my $wx (@tab) {
		my $ligne_taxo = $taxeuk{$wx};
		my @t= split(/\t/, $ligne_taxo);
		my $taxo = $t[1];
		print B "$a\t$wx\t$t[1]\n";
	
		$tax_hits{$a}{$i} = $t[1];
		$seq_hits{$a}{$i} = $tab[0];
		$i++;
	}
}	

close B;


my %seq_affil;

		if ( -e $NGS_id_Results."/panam_output/Phylogeny/Node_files/file_get_descendent_seq_$taxo_group"){
			open (W, $NGS_id_Results."/panam_output/Phylogeny/Node_files/file_get_descendent_seq_$taxo_group") ;

			my $i=0;
			my %hh; my %child; my %node_seq; my %boot;
			while (<W>) {
				if ($_=~ /Node :\s(.*?)\schild/) { 
					my $a = "Node_$i"; 
					my $boot_value = $1;
					my $child = ($_=~ tr/child//);
					$child{$a}= $child; 
					$child{$a}{'bootstrap'} = $boot_value;
	
					$boot{$a}= $boot_value; ## print "$a\t$boot_value\n";
					my @tab = split (/:/, $_);
					foreach my $e (@tab) { 
						my @t = split (/ /,$e);
						foreach my $p (@t) { 
							my $pp; 
							if (($p=~ /^\w/) and ($p !~ /\./)) { 
								chomp ($p);
								if (($p !~ /Node/) and ($p !~ /child/)){
									$pp = $p;
								}
								if (defined $pp ) {
									$node_seq{$pp} .="\t".$a;
									my @tax_seq = split(/\t/, $taxeuk{$pp});
									my $tax_seq = $tax_seq[0];
									if (defined ($tax_seq)) {
										chomp ($tax_seq);
										$hh{$a}{'tax'}.="\t".$tax_seq;
										$hh{$a}{'seq'}.="\t".$pp;	
	
									}	
								}
							}
						}
					}
				}
				$i++;
			}
	
	##############"NN	
			my %nn_node;
			foreach my $nn (keys %hh) {
				my @mmm = split (/\t/, $hh{$nn}{'tax'});
				if ($mmm[1] ne "") {
					$nn_node{$nn} = $mmm[1]
				}
			}
	
	###################"LCA
			my %taxo;
			foreach my $p (keys %hh) {	
				my @ttax=();
				my @tab = split (/\t/,$hh{$p}{'tax'}); 
				my $pp = $#tab;
				
			##### toutes les taxonomies reliées à un noeud sont stockées dans @t; 
		
				for (my $i=0; $i<=$pp; $i++) {
					if ($tab[$i] ne "") {
						push(@ttax, $tab[$i]);
					}
				}
				$taxo{$p}{'tax'} = \@ttax;
			}
	
			##########################

			my %taille_taxo;
			foreach my $bb (keys %taxo) {
				my $lmin_silva=100000; 
		
				my $mm = $#{$taxo{$bb}{'tax'}};
				my @r; my @rr; my $nnn ; my $nn;
				
				for(my $i=0; $i<=$mm; $i++) {
					my @r = split(/;/, ${$taxo{$bb}{'tax'}}[$i]);
					my $nn = $#r;
					if (($nn >=0) and ($nn < $lmin_silva)) {
						$lmin_silva = $nn;
					}
				}
					
				$taille_taxo{$bb}{'tax'} = $lmin_silva;
			}

			my %term; my %taxo_com; my $qq; my @tt; my %rang_taxo;
				
			foreach my $p (keys %taxo) {
				my $qs = $#{$taxo{$p}{'tax'}}; 
				for (my $oo =0; $oo<= $qs; $oo++) {
					##my @t = split (/\//, ${$taxo{$p}{'tax'}}[$oo]);
					my @t = split (/:/, ${$taxo{$p}{'tax'}}[$oo]);
					for (my $i = 0; $i <= $taille_taxo{$p}{'tax'}; $i++) {
						if (exists $t[$i]) {
# 							$term{$p}{'tax'}{$i}.="\t".$t[$i].";";
							$term{$p}{'tax'}{$i}.="\t".$t[$i];
						}
				
					}	
				}
			}
	
			my %tax_node;
			foreach my $tt (keys %term) {
				foreach my $cc (keys %{$term{$tt}}) {
					$tax_node{$tt}{$cc}=""; 
					my %temp=%{$term{$tt}{$cc}};
					my @liste_level=sort keys(%temp);
					foreach(@liste_level){
						my $ss=$_;
						my %rang_tax;		
						if (defined $term{$tt}{$cc}{$ss}) {
							my @ff = split(/\t/, $term{$tt}{$cc}{$ss});
							my $t_max = $#ff;
							foreach my $ll (@ff) {
								$rang_tax{$ll}++;
							} 
				
							foreach my $aa (keys %rang_tax) {
								if ($aa ne "") {
									if ($rang_tax{$aa} == $t_max) {
										#$tax_node{$tt}{$cc}.=$aa;
										$tax_node{$tt}{$cc}.=$aa.";";
									}
								}
							}
						}
					}
				}
			} 
	
			################ printing results 
			################ /!\ modif : panam_affiliation par grp taxo 28/3/2012
			open (R, ">".$NGS_id_Results."/".$panam_output."/PANAM_Affiliation_".$taxo_group.".txt");
			if (-e $NGS_id_Results."/".$panam_output."/Similarity_annotation/Sorted_Sequences/files_$taxo_group/file_cut") {
				open (G, $NGS_id_Results."/".$panam_output."/Similarity_annotation/Sorted_Sequences/files_$taxo_group/file_cut")  ;
				my %nom_seq;
				while (<G>) {	
					if ($_=~ /^>(.*?)\s+/) {
						my $aa = $1; 
						$nom_seq{$aa} = 1;
					}
				}

				foreach my $e (keys %nom_seq) { 
				##############"ajout 24/5	
					$seq_affil{$e}=1;
				###########fin ajout 24/5  
					if (exists $h{$e}) {
						my $min = 1000000;
						my $res_node;
						my @tab = split (/\t/, $node_seq{$e});
						foreach my $a (@tab) {
							if ($a ne "") {
								if (($child{$a} < $min) and defined ($nn_node{$a}) ) {
									$min = $child{$a}; 
									$res_node = $a; 
								}
							}
						}
	
						my @mmm = split (/\t/,$hh{$res_node}{'seq'});
						my $nearest = $mmm[1];
						my @ttax_nearest = split(/\t/,$taxeuk{$nearest});
						my $tax_nearest	= $ttax_nearest[0];
						
						my $tax_node;if ($tax_node{$res_node}{'tax'} eq "") { #print "res_node --- $res_node\n"; 
							$tax_node = "Other;"}	
						else { $tax_node = $tax_node{$res_node}{'tax'}	 }
				
						########## 

						my @ppp = split (/;/, $tax_nearest);
						my $ppp = $#ppp;
						if ($ppp == 1) {
							$tax_nearest .="unclassified ".$ppp[1].";"
						}
						
						my @ooo = split (/;/, $tax_node);
						my $ooo = $#ooo;
						if ($ooo == 1) { 
							$tax_node .="unclassified ".$ooo[1].";" ; 
						}
				
						########## /!\ 28/3/2012 print dans 2 fichiers : le final AFFIL et le provisoir $k
						print R ">Read: $e\tBootstrap:$boot{$res_node}\nNearest neighbor: $nearest\t$tax_nearest\nLowest node: $res_node\t$tax_node\n\n";
						
						#  Realise par concatenation etape 11
						#print AFFIL ">Read: $e\tBootstrap:$boot{$res_node}\nNearest neighbor: $nearest\t$tax_nearest\nLowest node: $res_node\t $tax_node\n\n";

					}
				}
			}
			close R;
		}

